package scu.maqiang.eigenvalue;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T3PE;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.T3Type;
import scu.maqiang.numeric.EigenSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class PEEigenValue2DT3 {

	public static void main(String[] args) {
		Mesh2T3 mesh = new Mesh2T3().square2D(100, 100, T3Type.Left);
		//mesh.scale(Math.PI);
		
		FES2T3PE fs = new FES2T3PE(mesh);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		SRMatrix M = new SRMatrix(fs.GetNdof());
		
		int N = 10;
		double[] D = new double[N];
		double[][] V = new double[N][fs.GetNdof()];
		double[] x = new double[fs.GetNdof()];
		Arrays.fill(x, 1.0);
		
//		double G1 = 1.0, ep1 = 0.5, e1 = 0.1;
//		fs.assemblePEStiff(new double[] {G1, ep1, e1}, BVPType.COMMON, A);
//		fs.assemblePEMass(new double[] {1.0}, BVPType.CONSIST_MATRIX, M);
		
		double[][] matD = new double[][] {{14.9e10, 9.8e10,  0.0,    0.0,     -5.4}, 
                                          {9.8e10,  14.3e10, 0.0,    0.0,     -5.4}, 
                                          {0.0,     0.0,     2.2e10, 11.7,    0.0}, 
                                          {0.0,     0.0,     -11.7,  8.0e-9,  0.0}, 
                                          {5.4,    5.4,        0.0,    0.0,     7.2e-9}};
        fs.assembleStiff(new double[][][] {matD}, BVPType.CONSTITUTE_MATRIX_COMMON, A);
		fs.assembleMass(new double[] {1.0e3}, BVPType.CONSIST_MASS, M);
		//M.sort();
		//System.out.println(M);
		//System.exit(0);
//		M.sort();
//		System.out.println(M);
//		System.exit(0);
		fs.applyBC_MBN(A, x, new int[] {0, 1}, new double[] {0.0, 0.0}, 1, 2, 3, 4);
	    fs.applyBC_MBN(A, x, new int[] {2}, new double[] {0.0}, 1, 2, 3, 4);
		//fs.applyBCZero(A, x, Direct.X, BCType.RCE, 1);
		//fs.applyBCZero(A, x, Direct.Y, BCType.RCE, 2);
		//fs.applyBCZero(A, x, Direct.X, BCType.RCE, 1, 2);
		//System.out.println(A);
		EigenSSolver solver = new EigenSSolver();
		solver.inversePowerTran_DSS_2(A, M, x, 100000.0, N, D, V);
		System.out.println(MVO.toString(D));
		//mesh.toTecplot("Eigenvalues1D.dat", V);
		String outputFilePreFix = "PEEigenValue2D";
		double[][] u = new double[3][mesh.getNv()];
		for(int i = 0; i < N; i++) {
			fs.extractComponent(V[i], u);
			mesh.toTecplot(outputFilePreFix + (i+1) + ".dat", u);
		}
	}

}
